using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Globalization;
using System.IO;

using Atomic.Libraries;
using Atomic.Libraries.Compression;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Mathematics.Statistics;
using Atomic.Phonons;
using Atomic.Structures;
using Atomic.Vasp.LinearResponse;
using Atomic.Thermodynamics.Phonons;

namespace Atomic.Interfaces.Yphon
{
	public static class YphonPhonon
	{
		public static PhononStateDensity CreateStateDensity(IYphonVaspLinearResponseResult result)
		{
			// Use unfiltered output by default.
			return CreateStateDensity(result, false, false, -1);
		}

		public static PhononStateDensity CreateStateDensity(IYphonVaspLinearResponseResult result, int mesh)
		{
			return CreateStateDensity(result, false, false, mesh);
		}

		public static PhononStateDensity CreateStateDensity(IYphonVaspLinearResponseResult result, int mesh, bool removeNegative, bool normalize)
		{
			if (mesh < 0)
			{
				throw new ArgumentOutOfRangeException();
			}

			return CreateStateDensity(result, removeNegative, normalize, mesh);
		}

		public static PhononStateDensity CreateStateDensity(IYphonVaspLinearResponseResult result, bool removeNegative, bool normalize)
		{
			return CreateStateDensity(result, removeNegative, normalize, -1);
		}

		private static PhononStateDensity CreateStateDensity(IYphonVaspLinearResponseResult result, bool removeNegative, bool normalize, int mesh)
		{
			bool debug = false;

			// Extract force constants from OUTCAR.
			IBinaryData superfij;
			YphonSuperfij.Generate(result.Contcar, result.Outcar, out superfij);

			using (TempDirectory temp = new TempDirectory())
			{
				// Yphon creates a bunch of temporary files. Don't mess with the original files.
				File.WriteAllBytes(temp.Directory.AddFile("superfij.out").FullName, superfij.GetData());

				Process p = new Process();
				p.StartInfo.FileName = "bash";
				p.StartInfo.WorkingDirectory = temp.Directory.FullName;
				p.StartInfo.UseShellExecute = false;
				p.StartInfo.RedirectStandardInput = true;
				if (!debug)
				{
					p.StartInfo.RedirectStandardOutput = true;
					p.StartInfo.RedirectStandardError = true;
				}

				p.Start();

				if (mesh == -1)
				{
					p.StandardInput.WriteLine("Yphon -bvec -thr2 0.01 < superfij.out");
				}
				else
				{
					p.StandardInput.WriteLine(string.Format("Yphon -nq {0} {0} {0} -bvec < superfij.out", mesh));
				}
				p.StandardInput.WriteLine("exit $?");
				p.StandardInput.Close();

				int atoms = 0;

				string line;
				while ((line = p.StandardOutput.ReadLine()) != null)
				{
					if (line.Contains("Atom(s) in the Primitive cell"))
					{
						atoms = int.Parse(line.Trim().Split(' ')[0]);
					}
				}

				if (atoms <= 0)
				{
					throw new Exception("Couldn't determine number of atoms in the primitive cell.");
				}

				p.WaitForExit();

				if (p.ExitCode != 0)
				{
					throw new Exception(string.Format("Failed running Yphon. Error message: {0}", p.StandardError.ReadToEnd()));
				}

				return ParseStateDensity(temp.Directory.AddFile("vdos.out"), atoms, removeNegative, normalize);
			}
		}

		public static YphonPhononDispersion CreateDispersion(VaspLinearResponseResult result, SpaceVector[] points)
		{
			return CreateDispersion(result, Transform(points));
		}

		public static YphonPhononDispersion CreateDispersion(VaspLinearResponseResult result, Tuple<SpaceVector, SpaceVector>[] segments)
		{
			throw new NotImplementedException();
			/*if (result.Phonons.Superfij == null)
			{
				throw new ArgumentException("Couldn't find output of vasp_fij.");
			}

			using (TempDirectory temp = new TempDirectory())
			{
				// Yphon creates a bunch of temporary files. Don't mess with the original files.
				result.Phonons.Superfij.Write(temp.Directory);

				// Branches along the edge of the irreducible Brillouin zone (must be specified; determined by the symmetries of the primitive cell).
				using (StreamWriter sw = new StreamWriter(temp.Directory.AddFile("dis").FullName))
				{
					foreach (Tuple<SpaceVector, SpaceVector> segment in segments)
					{
						sw.WriteLine(string.Format(CultureInfo.InvariantCulture, "{0:r} {1:r} {2:r} {3:r} {4:r} {5:r}", segment.Item1.X, segment.Item1.Y, segment.Item1.Z, segment.Item2.X, segment.Item2.Y, segment.Item2.Z));
					}
					sw.Close();
				}

				Process p = new Process();
				p.StartInfo.FileName = "bash";
				p.StartInfo.WorkingDirectory = temp.Directory.FullName;
				p.StartInfo.UseShellExecute = false;
				p.StartInfo.RedirectStandardInput = true;
				p.StartInfo.RedirectStandardOutput = true;
				p.StartInfo.RedirectStandardError = true;

				p.Start();
				p.StandardInput.WriteLine("Yphon -bvec -thr2 0.01 -pdis dis < superfij.out");
				//                        "Yphon -bvec -thr2 0.01 -plot  < superfij.out > aa; "
				//p.StandardInput.WriteLine("Yphon -symm fcc -bvec -thr2 0.01 -pdis dis < superfij.out");
				p.StandardInput.WriteLine("exit $?");
				p.StandardInput.Close();
				p.WaitForExit();

				if (p.ExitCode != 0)
				{
					throw new Exception(string.Format("Running Yphon failed. Error message: {0}", p.StandardError.ReadToEnd()));
				}

				int m = 3 * result.FinalStructure.PrimitiveStructure.Sites.Count;
				List<List<PhononDispersionPoint>> points = new List<List<PhononDispersionPoint>>();
				for (int i = 0; i < m; i++)
				{
					points.Add(new List<PhononDispersionPoint>());
				}

				double d = 0.0;
				SpaceVector q0 = null;
				foreach (string line in File.ReadLines(temp.Directory.AddFile("vdis.out").FullName))
				{
					if (line == "")
					{
						continue;
					}

					Vector a = AtomicEnvironment.ParseVector(line);
					SpaceVector q = new SpaceVector(a.GetVector(1, 3).ToArray());
					Vector f = a.GetVector(4, m);

					if (q0 != null)
					{
						d += SpaceVector.Norm(q - q0);
					}
					q0 = q;

					for (int i = 0; i < m; i++)
					{
						points[i].Add(new PhononDispersionPoint(new Position(result.FinalStructure.PrimitiveStructure.Lattice, q), d, f[i]));
					}
				}

				List<PhononDispersionBranch> bands = new List<PhononDispersionBranch>();
				for (int i = 0; i < m; i++)
				{
					bands.Add(new PhononDispersionBranch(points[i]));
				}

				File.AppendAllText(temp.Directory.AddFile("vdis.plt").FullName, "pause -1 \"Press return to continue...\"");

				// Store some original plotting data in addition to the parsed output. This allows
				// comparing the plots with the originally generated plot.
				return new YphonPhononDispersion(bands, new RawDataFile(temp.Directory.AddFile("vdis.out")), new RawDataFile(temp.Directory.AddFile("vdis.plt")), new RawDataFile(temp.Directory.AddFile("vline.dat")));
			}*/
		}

		public static YphonPhononDispersion CreateFaceCenteredCubicDispersion(VaspLinearResponseResult result)
		{
			SpaceVector[] points = new SpaceVector[] { new SpaceVector(1.0, 1.0, 1.0), new SpaceVector(0.5, 0.5, 1.0), new SpaceVector(0.25, 0.5, 0.75), new SpaceVector(0.5, 0.5, 1.0), new SpaceVector(0.0, 0.0, 0.0), new SpaceVector(0.5, 0.5, 0.5) };
			//string[] labels = new string[] { "G", "X", "W", "X'", "G", "L" };
			return CreateDispersion(result, points);
		}

		public static YphonPhononDispersion CreateHexagonalDispersion(VaspLinearResponseResult result)
		{
			throw new NotImplementedException();
		}

		private static PhononStateDensity ParseStateDensity(FileInfo file, int primitiveStructureSize, bool removeNegative, bool normalize)
		{
			List<PhononStateDensityPoint> points = new List<PhononStateDensityPoint>();
			MeanStatistics delta = new MeanStatistics();

			foreach (string line0 in File.ReadLines(file.FullName))
			{
				string line = line0.Trim();
				while (line.Contains("  "))
				{
					line = line.Replace("  ", " ");
				}

				if (line != "")
				{
					string[] cols = line.Split(' ');
					double frequency = double.Parse(cols[0], CultureInfo.InvariantCulture);
					double density = double.Parse(cols[1], CultureInfo.InvariantCulture);

					if (frequency <= 0.0 && removeNegative)
					{
						continue;
					}

					if (points.Count > 0)
					{
						delta.Add(frequency - points[points.Count - 1].Frequency);
					}
					points.Add(new PhononStateDensityPoint(frequency, density));
				}
			}

			double delta0 = delta.Mean;

			if (normalize)
			{
				double integral = 0.0;
				for (int i = 0; i < points.Count; i++)
				{
					integral += delta0 * Math.Max(0.0, points[i].Density);
				}

				for (int i = 0; i < points.Count; i++)
				{
					points[i] = new PhononStateDensityPoint(points[i].Frequency, Math.Max(0.0, points[i].Density) * 3 * primitiveStructureSize / integral);
				}
			}

			return new PhononStateDensity(points, 3 * primitiveStructureSize);
		}

		private static Tuple<SpaceVector, SpaceVector>[] Transform(SpaceVector[] points)
		{
			if (points.Length < 2)
			{
				throw new ArgumentException("Invalid number of phonon branch points");
			}

			List<Tuple<SpaceVector, SpaceVector>> branches = new List<Tuple<SpaceVector, SpaceVector>>();
			for (int i = 1; i < points.Length; i++)
			{
				branches.Add(Tuple.Create<SpaceVector, SpaceVector>(points[i - 1], points[i]));
			}

			return branches.ToArray();
		}

		/*

		public PhononDispersionCurve CreateFaceCenteredCubicDispersionCurve()
		{
			Vector3[] points = new Vector3[] { new Vector3(1.0, 1.0, 1.0), new Vector3(0.5, 0.5, 1.0), new Vector3(0.25, 0.5, 0.75), new Vector3(0.5, 0.5, 1.0), new Vector3(0.0, 0.0, 0.0), new Vector3(0.5, 0.5, 0.5) };
			string[] labels = new string[] { "G", "X", "W", "X'", "G", "L" };
			return PhononDispersionCurve.Generate(Superfij, points, labels);
		}

		public PhononDispersionCurve CreateHexagonalDispersionCurve()
		{
			Vector3[] points = new Vector3[] { new Vector3(0.0, 0.0, 0.0), new Vector3(0.0, 0.0, 0.5), new Vector3(0.5, 0.0, 0.0), new Vector3(0.5, 0.0, 0.5), new Vector3(1.0 / 3, 1.0 / 3, 0.0), new Vector3(0.0, 0.0, 0.0) };
			string[] labels = new string[] { "K0", "K1", "K2", "K3", "K4", "K5" };
			return PhononDispersionCurve.Generate(Superfij, points, labels);
		}*/
	}
}
